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Solar heating of acometary surface provides the energy necessary to sustain gaseous 
activity, through which dust is removed’. In this dynamical environment, both the 
coma** and the nucleus”® evolve during the orbit, changing their physical and 
compositional properties. The environment around an active nucleus is populated by 
dust grains with complex and variegated shapes’, lifted and diffused by gases freed 
from the sublimation of surface ices®”. The visible colour of dust particles is highly 
variable: carbonaceous organic material-rich grains” appear red while magnesium 
silicate-rich’’” and water-ice-rich’” grains appear blue, with some dependence on 
grain size distribution, viewing geometry, activity level and comet family type. We 
know that local colour changes are associated with grain size variations, such as in the 
bluer jets made of submicrometre grains on comet Hale-Bopp” or inthe fragmented 
grains in the coma” of C/1999 S4 (LINEAR). Apart from grain size, composition also 
influences the coma’s colour response, because transparent volatiles canintroducea 
substantial blueing in scattered light, as observed in the dust particles ejected after 
the collision of the Deep Impact probe with comet 9P/Tempel 1”. Here we report 
observations of two opposite seasonal colour cycles inthe coma and on the surface of 
comet 67P/Churyumov-Gerasimenko through its perihelion passage’®. Spectral 
analysis indicates an enrichment of submicrometre grains made of organic material 
and amorphous carbon inthe coma, causing reddening during the passage. At the 
same time, the progressive removal of dust from the nucleus causes the exposure of 
more pristine and bluish icy layers on the surface. Far from the Sun, we find that the 
abundance of water ice on the nucleus is reduced owing to redeposition of dust and 
dehydration of the surface layer while the coma becomes less red. 


Understanding how comets work and evolve is one of the most compel- 
ling questions to which the Rosetta mission’’ has been trying to find an 
answer. Unlike past fly-by missions to comets, Rosetta was designed 
to enter the same orbit and accompany comet 67P/Churyumov-Ger- 
asimenko (hereafter 67P) through its perihelion passage, giving us 
the unique opportunity to follow the colour changes developing on 
a comet nucleus and coma in its active phase through the inner Solar 
System. So far, several studies have investigated how dust and gas 
(H,O, CO.) production’ are correlated during one cometary rotation, 
showing that the water-ice sublimation on the surface of 67P is the 
driving mechanism for dust ejection in the coma7. In fact, the dust 
emission flux appears lower above high-cohesion consolidated terrains 
and reaches amaximum when the subsolar direction is aligned above 
volatile-rich areas. As the illumination conditions continuously change 
on 67P, rather than focusing on coma and nucleus colour variations 
occurring during a few diurnal rotations, we analysed the entire dataset 


returned by the Visible, Infrared and Thermal Imaging Spectrometer 
(VIRTIS)”’ on Rosetta. Our method allows exploration of the comet’s 
colour evolution starting fromJanuary 2015 (inbound orbit, heliocen- 
tric distance 2.55 AU), encompassing perihelion passage (August 2015, 
1.24 AU) until May 2016 (outbound orbit, 2.92 AU). Spectral indicators in 
the visible range, including the integrated radiance (/), the wavelength 
of the radiance peak (A,,,,,) and the spectral slopes (in the 0.4-0.5 and 
0.5-0.8 um ranges), are synergistically used to model the composition 
and grain size distribution of the dust particles in the coma as a func- 
tion of heliocentric distance. The spectral indicators are derived for 
each spectrum acquired by VIRTISin an annulus encircling the nucleus 
and comprising all the pixels inthe coma located at a tangent altitude 
between 1.0 and 2.5 km above the nucleus’ limb. At the same time, we 
monitor the colour changes occurring on 12 test areas on the nucleus 
through the 0.5-0.8 pm spectral slope**. The description of the dataset, 
including an example of VIRTIS images and spectra after calibration 
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reduced dataset, 


complete annulus coverage (red 
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Fig. 1| Time series of the spectral 


properties of 67P coma dust. 


a, Integrated radiance (/): full 
points; see discussion in Methods). 


b, The 0.4—0.5 um spectral slope. 
c, Visible radiance wavelength peak 


Anax: d, The 0.5-0.8 pm spectral 


dataset, partial annulus coverage 
slope. Throughout the paper, 


spectral slopes aregiveninpm” 


(black points); 
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corresponding to 10%/100 nm colour 


slope. Rosetta’s MTPintervals are 
marked by vertical lines between 
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and heliocentric distance (D,) of MTP 
periods are listed in Extended Data 


Table 1. Perihelion occurs during 
MTPO19. Dates are displayed as year- 


MTPO12 and MTPO28. Time intervals 
month-day. 
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regionname and position (longitude, lon, and latitude, lat) of each area are 
reported above each plot. The maximum error associated with the slope 


measurements is of the order of +O.1jum™. 


Fig. 2| Time evolution of 67P nucleus colour measured through the 0.5- 


0.8 xm spectral slope above 12 control areas. Perihelion passage occurs 


during MTPO19 and is marked by the vertical dashed line. The morphological 
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Fig. 3 | Simulations of the /,,,, as a function of the scattering angle for 
spherical particles of different composition. a, Amorphous carbon. 

b, Magnesium silicate. c, Organic ice tholin (A,,,, > 0.65 pm for 0.10-0.99-um- 
radius grains). d, Troilite. e, Water ice. f, 67P observations by VIRTIS. Median 


improvement, spectral indicator processing and data modelling, are 
given in Methods. 

From a synergistic time analysis of the coma and of the nucleus sur- 
face spectral properties, two opposite trends with heliocentric distance 
are observed. Coma integrated radiance increases (See/time series in 
Fig. 1a) and colour becomes redder when the comet approaches peri- 
helion (see time series of /,,,,, and the 0.4—0.5 and 0.5-0.8 pum spectral 
slopes in Fig. lb-d). Conversely, on12 control areas of the nucleus, for 
which we have continuous time coverage during the entire Rosetta 
mission, we observe a systematic blueing of the surface (see 0.5—-0.8 
Lum spectral slope time series in Fig. 2), with a reduction of the spectral 
slope up to 50% measured between the passage of 67P through the 
frost line (heliocentric distance 2.7 AU, occurring in October 2014 dur- 
ing inbound orbit and in April 2016 on the outbound) and perihelion 
(1.25 AU, August 2015). During the pre-perihelion period, we interpret 
the coma colour change as being the consequence of the progressive 
loss of the ice fraction in the dust grains ejected from the nucleus”, 
which makes them redder. 

When the comet approaches the Sun, the dust grains lifted from the 
nucleus’ southern hemisphere—the one illuminated at that time—are 
warmer and become more dehydrated. Simulations based on Mie scat- 
tering theory applied to spherical grains of radius between 0.land50 pm 
and five different compositions (see Methods for a discussion about the 
scattering model limitations induced by spherical shaped grains) are used 
tomodel VIRTIS spectral indicators, for example, A,,,, and the 0.5—0.8 pm 
spectral slope. Simulation results (Figs. 3, 4) are compatible with the 
presence of submicrometre grains made of organic matter and carbon 
during perihelion passage, which could be the result of the fragmentation 
of the mineral and organic matrix embedded in larger grains. During the 
outbound orbit, the simulations show that the organic material contin- 
ues to be present in the coma together with an increasing percentage 
of different grains (water ice or magnesium silicate, both of which are 
not fully compatible with the spectral indicator trends (Extended Data 
Table3)), contributing to the progressive blueing observed in the data. In 
this orbital phase, the Sun shines again on thenorthern hemisphere, from 
where a greater flux of surficial water-ice-rich grains” is lifted inthe coma. 
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phases. Standard deviation computed on VIRTIS data are of the order of 
0.03 um. 


Conversely, the nucleus’ surface becomes progressively bluer during 
the inbound orbit because the greater solar input causes an increase in 
the gaseous activity. This triggers dust ejection with consequent ero- 
sion of the surface: on average, a0.5-m-thick erosion layer is lost during 
each orbit”. As a result of this process, more pristine subsurface layers 
enriched in water ice are exposed on the surface, which turns bluer?*”. 
We have evidence*®, in fact, that water-ice-rich layers are immediately 
available under the dust layer, as shown in areas of recent disruption 
of the surface” and as demonstrated by many thermal models applied 
to cometary nuclei?***. Between 20% (ref. *) and 50% (ref. 7°) of the 
dust grains ejected around perihelion from the south pole—the region 
that receives the maximum insolation at that time—fall back prefer- 
entially on the northern hemisphere. The fall-back flux is dominated 
by decimetre-size dust aggregates in which a small fraction of water 
ice can be maintained”. After perihelion, as soon as activity begins 
to settle, the progressive accumulation of dust on the surface covers 
the exposed pristine layers again, making the nucleus surface redder 
again. This process explains the colour cycle observed by VIRTIS on 
the nucleus (Fig. 2). 

The colour of the 67P coma particles can be modelled with different 
grain populations made of water ice before perihelion, organic mate- 
rial and carbon at perihelion, and water ice or possibly magnesium 
silicate after perihelion. Similar composition endmembers are com- 
patible with the findings reported by previous Rosetta studies: during 
outbursts, the release of submicrometre-size water-ice grains mixed 
with hundred-micrometre-size refractory grains has been observed”. 
Atthe same time, a broad emission feature between 3.4 and 3.6 um has 
been detected by VIRTIS”®, which is compatible with the sublimation 
of organic materials caused by the high temperatures of the grains. 
The colour reduction observed when the comet was far from the Sun 
resembles the spectral properties of water-ice grains, whose existence 
in the lower coma has been confirmed by other investigations: apart 
from the presence of submicrometre-size water-ice grains ejected 
during outbursts”, thermal models and Optical, Spectroscopic, and 
Infrared Remote Imaging System (OSIRIS) camera images”’ show that 
at 1.53 AU heliocentric distance, the total sublimation of ice for grains of 
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Fig. 4| Simulations of the 0.5-0.8 um spectral slope as a function of the 
scattering angle for spherical particles of different composition. a, 
Amorphous carbon.b, Troilite.c, Magnesium silicate. d, Water ice.e, Organics. 
f, Median spectral slope derived from VIRTIS observations for pre-perihelion 
(MTPO15 and MTPO16), perihelion (MTPO20 and MTPO21) and post-perihelion 
(MTPO27 and MTPO28) mission phases. The standard deviations of the VIRTIS 
dataare<0.6um". 


radius between 5 and 500 pm made of ice-dust layers or mixtures occurs 
well beyond the 2.5 km distance from the 67P nucleus considered here. 
This implies that water-ice grains can populate the region of the coma 
analysed in this study. Finally, a photometric investigation of the visible 
colours of more than1,000 individual grains resolved on OSIRIS colour 
images of the 67P coma”’ indicates the presence of three composition 
classes: organic material grains (steep colour slope), mixtures of silicate 
and organic material (intermediate slope), and water ice (flat slope). 
The presence of carbon, silicates and magnesium in 67P dust grains has 
been assessed by the Cometary Secondary Ion Mass Analyser (COSIMA) 
instrument on Rosetta’’. 

The fact that the coma’s and the surface’s colour spectral indica- 
tors return to the same values at the beginning and at the end of their 
time series (Fig. 1), for example, at heliocentric distances greater than 
2.7 AU, is evidence of the establishment of an orbital water-ice cycle 
driven by the solar heating. Far from the Sun, with the settling down 
of the gaseous activity’, a tenuous (/= O in Fig. 1a) coma surrounds a 
red-coloured dusty nucleus. 
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Methods 


VIRTIS dataset 

While the scientific objectives of the Rosetta mission”, including the 
study of the dust in the coma, were well-defined before the spacecraft 
reached 67P, the detailed planning of the scientific payload opera- 
tions has been a complex and continuously evolving process due to 
the uncertainties of the cometary environment in which the mission 
had to operate. A detailed overview of the Rosetta science planning 
process and challenges is given in ref. *. In this context, VIRTIS coma 
observations were performed from very different viewing and illumi- 
nation geometries, spacecraft ranges from the nucleus, heliocentric 
distances and cometary activity levels, resulting in a dataset showing 
avery variable spatial resolution and signal-to-noise ratio. 

While VIRTIS-led targeted observations of the coma have been 
routinely executed, a large fraction of the VIRTIS dataset consisted of 
ride-along and serendipitous acquisitions that were acquired while 
the spacecraft attitude and pointing were controlled by other instru- 
ments. For this reason, to maximize the retrieval of information about 
the dust particles properties, a statistical approach has been applied to 
the vast VIRTIS dataset with the scope to exploit itin its completeness. 

We restrict our analysis to only the VIRTIS-M (imaging channel) vis- 
ible data covering the 0.25-1.0 um spectral range with 432 bands anda 
spectral sampling of less than 2 nm per band. A detailed description of 
the VIRTIS experiment is provided by ref. ”. In contrast to the VIRTIS-M 
infrared channel operating within the 1-5 pm spectral range, which 
ceased operation in May 2015 due to a permanent failure of the cryo- 
cooler cooling of the infrared channel detector, the visible channel has 
operated flawlessly for the entire duration of the mission. 

For this study, we processed coma data collected during Rosetta’s 
MTPO12-MTPO28 periods. The time intervals of each medium-term 
phase (MTP) are listed in Extended Data Table 1. The dataset consists of 
more than 4,500 hyperspectral cubes of the coma of 67P acquired froma 
range of distances from the nucleus between 27.8 and 1,462.2 km, result- 
ing ina spatial resolution between 7 and 365 m per pixel, respectively. 

Each VIRTIS observation has been calibrated in spectral radiance” 
and specific geometry parameters have been calculated for each pixel 
falling on the coma region and the nucleus. Further details about the 
VIRTIS data mapping methods are given in ref. >. The nucleus solar 
phase and the spacecraft position and distance in the Cartesian frame 
centred on the nucleus have been calculated using SPICE kernels”. 

The spectral parameters used in our analysis to study the coma emis- 
sion are calculated in an annulus defined by a tangent altitude running 
between q,,,, =1kmand d,,,, =2.5 km around the nucleus. This criterion 
is adopted independently from the nucleus-spacecraft distance and the 
solar phase angle at the time of the observation. The selection of pixels 
placed at distances larger than1 km allows (1) reduction of the contribu- 
tion of the instrumental stray light coming from the illuminated part of 
the nucleus and (2) maximization of the coverage across the dataset. 

An example of a VIRTIS-M visible image of the coma showing the 
annulus definition is given in Extended Data Fig. 1 (top panels). This par- 
ticular image refers to observation MTPO19 - STP068 - V1_ 00397724286 
acquired on 9 August 2015 froma distance of 304 km, resulting in a spa- 
tial resolution of 75 m per pixel with an integration time of 16s per line. 

Sucha viewing geometry is not very common because the complete- 
ness of the annulus over the 360° azimuth (up to distance q,,,,) is limited 
by the instrumental field of view (FOV = 3.6°), by the relative distance 
(d) between the Rosetta spacecraft and the comet’s surface, and by the 
spacecraft off-nadir pointing direction (a). Only when the following 
condition is verified, VIRTIS-M can acquire the full annulus: 
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wherer=1.74 kmis the average radius of the nucleus”. A similar condi- 
tionis verified in only 248 observations (marked as red points in Fig. 1a), 


corresponding to about 5% of the dataset: as aconsequence, integrated 
radiance and maximum emission wavelength values could be biased 
by the incompleteness of the annulus on the remaining observations. 
The average integrated radiance on incomplete annulus observations 
(black points) appears, in general, higher than on the ones where the 
annulusis entirely acquired (red points). This happens because VIRTIS 
is preferentially observing the coma towards the solar direction to 
maximize the signal-to-noise ratio. Apart from the above limitations, 
which introduce some scattering in the data points, the value of the 
integrated radiance is mainly driven by the number of scatterers (dust 
grains) along the line of sight, rather than by their composition and 
grain size distribution (which the spectral slopes and the spectral posi- 
tion of the radiance peak depend on). 

The overall coma brightness and spectral properties measured by 
VIRTIS are driven by light scattering on dust particles and are not altered 
by emissions from gas molecules in the coma. The instrument sensitiv- 
ity at visible wavelengths, in fact, is not enough to measure gaseous 
emissions from CN, C, and C,. 


Calibration update using stellar observations 
With the aim of improving the VIRTIS-M visible response in the 0.25- 
0.45 um spectral range, where the standard pipeline suffers large uncer- 
tainties due to stray light introduced by the on-ground calibration setup, 
we applied a radiometric correction based on stellar observations. The 
average signal derived from Vega (a Lyrae) observations listed in Extended 
Data Table 2 is used as a reference to correct the standard responsivity. 
Four observations of Vega acquired with the longest possible integra- 
tion time, 50s, were used to retrieve the raw signal above the noise level 
(about +2 datanumbers (DN)) estimated on the nearby deep-sky pixels. 
Before this step, each observation was dark-level-subtracted and 
despiked (Extended Data Fig. 2 (top-left panel), showing the signal and 
sky level for one of the four Vega observations). Owing to the instrumen- 
tal point spread function and spectral tilt, both described in detail in 
ref.*>, the Vega spectral signal is distributed across 18 spatial pixels and 
is averaged above these. The four observations were processed in the 
same way before being further averaged to improve the signal-to-noise 
ratio. The resulting average Vega signal is shown in Extended Data Fig. 2 
(top-right panel). The fractional deviation of the single spectra fromthe 
average is about 0.05. The Vega flux in ref. *° is convolved with the VIRTIS 
Gaussian response (full-width at half-maximum, 2.3 nm) and then bya 
4-nm-wide boxcar. The spectra were then placed on an absolute wave- 
length scale using the centre of Balmer and Paschen intense absorption 
lines as references and assuming a constant dispersion. The spectral 
dispersion (in nm) applied is A(b) = 232.9 + 1.88515 x b, with b ranging 
from Oto 431, corresponding to the VIRTIS-M number of spectral bands. 
The resulting relative responsivity derived as the ratio between VIRTIS 
signal and Vega flux is shown in Extended Data Fig. 2 (bottom-left panel). 
To further remove residual high-frequency noise due to the low Vega 
signal, anine-band running boxcar filter was applied. The response 
curves were normalized at the value of 58.75 at A = 0.635 tm to allow 
a direct comparison with the standard pipeline responsivity shown 
in Extended Data Fig. 2 (bottom-right panel). Within the limits of the 
low Vega signal in the blue spectral range (less than 40 DN at 0.55 pm), 
this method allows the retrieval of a better responsivity at shorter 
wavelengths with respect to the standard pipeline. As shown by the 
ratio plot (blue curve), the standard pipeline overestimates the target 
radiance for wavelengths A < 0.4 pm. In the remaining spectral range, 
the differences are much smaller. As the wavelength of the dust emis- 
sion spectral radiance peak occurs between 0.45 and 0.55 pm, such an 
instrument-response correction allows improvement of the retrieval of 
the left wing of the radiance peak and the spectral slope determinations. 


Coma spectral indicators 
After averaging the coma signal within the annulus previously defined 
and applying the updated responsivity function, we calculated four 
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spectral indicators for each observation. The indicators are as fol- 
lows. (1) The integrated radiance (/) across the visible spectral range, 
defined as: 


lum 


Ent S orsum RAGA 
[= Oy 


(1) 


where R(n,A) is the spectral radiance measured on the nth pixel of the 
annulus at wavelength A and Nis the total number of pixels within the 
annulus area. (2) The wavelength of maximum emission of the radiance 
(A na.) Measured ona fourth-degree fit on the average spectral radiance 
calculated on the annulus. 

After converting spectral radiance to irradiance/solar flux (//F), 
we calculated the two spectral slopes following ref. °. These are: (3) 
0.4—0.5 um spectral slope on//F; and (4) 0.5—0.8 um spectral slope on 
I/F. The two spectral slopes are calculated as the angular coefficient 
of the best linear fit in the ranges 0.4-0.5 and 0.5-0.8 pm after having 
normalized //F at 0.5 um. The determination of the spectral indicators 
follows the scheme shown in Extended Data Fig. 1 (bottom panels). 

Spectral indicators depend on illumination geometries and the 
spacecraft—nucleus distance: as a consequence of the demanding mis- 
sion scenario, complex orbits, including terminator, icosahedron, 
flybys and far excursions, have been flown”, resulting in a wide range 
of distances and solar phases. The correlation between the integrated 
radiance /and distance is shown in Extended Data Fig. 3, where VIRTIS-M 
measurements are co-located with Rosetta’s position and solar phase 
angle with respect to the comet nucleus along the flight trajectory. 


Temporal evolution of the coma colour 

The integrated radiance (/) time series (Fig. la) shows an increase in 
the coma’s brightness starting from January 2015 (MTPO12, heliocen- 
tric distance 2.55 AU) to perihelion passage (MTPO19, 1.24 AU, August 
2015). Owing to the characteristics of Rosetta’s orbit around the nucleus 
and the different observation strategies implemented during the mis- 
sion, the completeness of the annulus region on the VIRTIS dataset is 
reached on only a limited number of observations, marked as red points 
in Fig. la. For the majority of the remaining observations (black points), 
partial coverage of the annulus is achieved. Notwithstanding the rela- 
tive scatters of the data introduced by the variability of the observ- 
ing conditions, some clear trends are identified. A local maximum, 
25</<30Wm’°sr", is measured about one month after perihelion. 
Two additional intensity peaks are measured at the end of August and 
September 2015, when a maximum intensity of about 50 Wm” sr‘ is 
reached. These peaks are caused by energetic cometary activity—two 
outbursts observed by VIRTIS on 13-14 September 2015 are reported 
by refs. 7°°’—but appear also to be strongly correlated with the peculiar 
viewing geometry and spacecraft position with respect to the nucleus: 
in fact, at the time of these observations, Rosetta was performing far- 
nucleus excursions, orbiting at greater distances than usual from the 
nucleus, up to 450 and 1,462 km, respectively. As aconsequence of the 
increased distance from the nucleus, VIRTIS-M has recorded higher 
radiances boosted by the larger column density of dust particles along 
the line of sight and by the maximum activity occurring immediately 
after perihelion. After this period, the integrated radiance drops rap- 
idly from MTPO21 to MTPO24 where an average value of 5W m” sr‘ is 
reached. The maximum integrated radiance measured in MTPO20 and 
MTPO21 is about a factor of 10 more intense than the average value of 
5 Wm” sr’ measured when the comet was far from the Sun (greater 
than 2 AU) on both inbound and outbound legs of the orbit. With the 
exclusion of the two peaks, a similar asymmetric, or cusp-like trend 
with amaximum occurring about one month after perihelion, has been 
measured by other Rosetta’s instruments for the water production 
rate, which is associated with the ejection of the dust®. Apart fromthe 
greater column density, we have clues that a larger particle albedo 
could contribute to the increase of the integrated radiance: according 


to ref.*, the albedo of the dust particles increases up to 20% near peri- 
helion for observations acquired at a 90° solar phase angle. Similar 
albedo changes could be the result of a different composition of the 
dominant grain population in the coma along the orbit, as discussed 
in the main text. 

The time series of the wavelength of maximum emission (A,,,,) 
(Fig. lc) shows a cusp-like trend starting from low values dispersed 
around J,,,,, = 0.45 um on January 2015 (MTPO12) on the inbound leg 
of 67P trajectory to/,,,, > 0.5 pm during perihelion passage (MTPO19). 
Moving along the outbound leg of the orbit, the wavelength of maxi- 
mum emission shifts progressively towards shorter wavelengths, 
reaching again A,,,, = 0.45 pm in mid-May 2016 (MTPO28). Whereas 
Ana, and integrated radiance time series show a similar time trend, 
spectral slopes time series (Fig. 1b, d) are characterized by a different 
evolution: at the two extremes of the time series, both 0.4—0.5 pm and 
0.5-0.8 um slopes are almost neutral whereas at perihelion they reach 
values dispersed at 5 and 3 pm’, respectively. The visible colours of 
the coma show asystematic reddening before perihelion passage and 
a progressive blueing after it. The concurrent variations occurring on 
all spectral indicator time series point to a change in the dust grain 
composition and grain size distribution with the heliocentric distance. 


Temporal evolution of the nucleus colour 

When observed at a global scale, 67P surface dust appears very dark 
(geometric visible albedo 6.2%; ref. °°) and dehydrated, with an average 
water-ice fraction of about 1% dispersed in submicrometre grains”. 
The residual water-ice content, if present, on the dust grains ejected 
fromthe surface undergoes rapid sublimation once the grainis lifted”. 
So far, infrared spectroscopic observations have confirmed the pres- 
ence of only two ices on the surface of 67P nucleus: water ice” and 
carbon dioxide”. 

Thelow thermal capacity of the dust (700Jkg“K™") andthe high solar 
flux can rapidly increase the temperature of the grains causing the 
sublimation of the volatiles. A colour temperature as high as 630 K is 
measured on dust grains during outbursts”’ and in the 260-320 K range 
on quiescent coma*. Furthermore, the grains observed by Rosetta’s dust 
instruments appear dehydrated*°”’. This probably happens because 
the instruments aboard the spacecraft are not kept at cryogenic tem- 
peratures, allowing the dispersion of the volatile fraction of the grains 
after their collection. Besides dehydration, grain size distribution also 
plays a role in the changes of the visible colours (see ‘Modelling dust 
grain composition and size distribution in the coma’). 

During the pre-perihelion phase, the nucleus’ surface shows a pro- 
gressive blueing until one month after perihelion when the minimum 
slope is reached in MTPO2O0 (Fig. 2). In this phase, the gas activity 
removes the dust from the surface with greater efficiency, resulting in 
the exposure of more pristine subsurface material enriched inices*°”. 
This trend is interrupted after the perihelion passage when a progres- 
sive surface reddening is observed. Along the outbound trajectory, 
the gaseous activity progressively settles, resulting in the accumula- 
tion of dehydrated (red) dust on the surface. This appears to bea gen- 
eral mechanism, occurring on the majority of the nucleus surface, as 
shown in Fig. 2, where the average 0.5—0.8 pm slope is shown for twelve 
30 m x 30-m-wide control areas at one-month time resolution. These 
control areas, selected between latitude +45.5° and —34.5°, are the 
ones that have been observed in daylight by VIRTIS during the entire 
duration of the mission, allowing us to follow the seasonal colour cycle 
developing on the nucleus. 


Dust-scattering simulations 

Thesimulations are computed following the Bohren—Huffman Mie scat- 
tering code to calculate scattering and absorption by ahomogenous 
isotropic sphere. The method relies on the following assumptions. 
(1) Single-scattering approximation, for example, the solar photons 
interact with only one particle before reaching the observer. Such a 


hypothesis is reasonable for the annulus region considered in this 
study, where the optical depth is small. (2) Polydisperse spherical grain 
size distributions with radii of 0.10-0.99 um (90 bins, 0.01 pm per 
bin), 1.0—9.9 pm (90 bins, 0.1 pm per bin) and 10-50 um (40 bins, 1 um 
per bin). (3) Five different grain compositions are simulated: water 
ice (optical constants from ref. *), organic material”®, troilite*, amor- 
phous carbon” and magnesium silicate*®. These endmembersare rep- 
resentative of the cometary nuclei composition, which is made up of 
a macromolecular assemblage in which various organic components 
(both aromatic and aliphatic), minerals (including silicates, iron sul- 
phides like pyrrhotite and/or troilite, and possibly ammoniated salt) 
and ices are mixed together™. (4) Spectral simulations are performed 
at 27 wavelengths within the 0.35-1.0 pm spectral range. (5) Scattering 
angle 0 (=180° — g, where gis the solar phase angle) is computed ata 
step of 1° from 55° to 155° and then averaged on 4° bins to match the 
VIRTIS FOV (3.662). 
Following ref. °°, the single scattering intensity is given by: 


KA) = SU) Rea(p)0q(a) PS (2) 


where S(A) is the solar flux” scaled at the comet’s heliocentric distance, 
No is the dust column density calculated along the line of sight from 
an observer placed at distance p from the nucleus, ois the grain geo- 
metrical cross-section, g(A) is the scattering efficiency of a single par- 
ticle and p(g, A) is the scattering phase function calculated at phase 
g. Rather than modelling the absolute value of the intensity, we limit 
here our analysis to the changes occurring in the spectral response of 
the average radiance emitted from the dust. In particular, we focus on 
the changes of the wavelength of the maximum emission (Fig. 1c) and 
0.5-0.8 um spectral slope (Fig. 1d) observed at different heliocentric 
distances. To achieve this goal, we calculated the spectral radiance fac- 
tor, (A) expressed in arbitrary units, derived from Eq. (2) by removing 
the factors depending on the dust grain distribution and properties: 


(A, g) « s(ayqiay PB” (3) 


The phase functions show scattering characteristics depending on 
composition and grain sizes: as a general rule, transparent grains, like 
water ice, are characterized by avery intense and collimated forward- 
scattering peak and complex resonance peaks at intermediate phase 
angles. Opaque materials, like troilite, show an almost isotropic scat- 
tering function. Semi-transparent particles, like ice tholins, have an 
intermediate behaviour witha prevalence of backscattering response. 
Grain size plays a fundamental role inthe scattering mechanism: inthe 
Rayleigh regime (particle size much smaller than the wavelength), light 
is equally forward- and backscattered while in the geometric regime 
(particle size much larger than the wavelength), backscattering domi- 
nates. 

Fixing composition and grain size, the corresponding phase func- 
tions at the 27 visible wavelengths allow the derivation of the spectral 
radiance according to Eq. (3). 

While p(g, A) has a smooth response in the Rayleigh and in the geo- 
metric scattering regimes, the simulations, in general, show more fluc- 
tuations occurring for certain combinations of composition, grain size 
and @ angle when the grain sizes are comparable to the wavelength, 
for example, at about 1 pm. To mitigate this effect, we computed 
simulations for 25° < g< 125°, corresponding to scattering angles of 
55° < @< 155°, at steps of 1° and then averaged the angular response 
of p(g, A) on 4° bins, similar to the VIRTIS FOV. Similarly, the p(g) was 
computed for single grain radius (0.10-0.99 ppm (90 bins, 0.01 pm per 
bin), 1.0—9.9 pm (90 bins, 0.1 pm per bin) and 10-50 um (40 bins, lum 
per bin)) and then averaged for each of the three families. Following this 
method, the /(/,g) (from Eq. (3)) was derived for each composition and 


size distribution. The A,,,,, and, after conversion in //F, the 0.5-0.8 um 
spectral slope were derived from the simulated /(A,g). The correspond- 
ing values are shown in Figs. 3, 4, respectively. 


Modelling dust grain composition and size distribution inthe 
coma 

To disentangle the viewing and illumination geometry effects from 
particles’ physical properties and to provide a more quantitative assess- 
ment of the changes observed in the A,,,,, time series, the light scattering 
on grains of different composition and radius was simulated. Optical 
and in situ measurements of the dust grain distribution in the 67P coma 
by Rosetta’s instruments have revealed that the particles lie mainly 
within the 0.1 pm-1 mm diameter range”. For particles smaller than 
1mm, the dust size distribution follows a power-law distribution witha 
coefficient varying between —2 for heliocentric distances beyond 2 AU 
and —3.7 at perihelion. Conversely, large grains (diameter greater than 
1mm) showa muchsteeper distribution with a coefficient of —4 (ref. >’). 
As small grains are the predominant population, in our analysis we 
preferentially model the spectral behaviour of small (radius less than 
50 um), compact grains. A similar grain size range has been also used 
by refs. *°°” to model dust optical properties at infrared wavelengths 
and by ref. *° to study thermal evolution and depletion of ice in dust 
aggregates along jets. The VIRTIS data are interpreted through the Mie 
theory, which can simulate the scattering properties of spherical grains 
of radius Rand homogeneous composition having x= (2TtR/A) < 1,000. 
Such a limit corresponds to a maximum grain radius R=50 pminthe 
visible spectral range (wavelength 0.35 <A < 1.0 pm). To simulate the 
spectral response of larger grains, including fluffy aggregates”, it is 
necessary to use a geometric optics approximation”, which is beyond 
the scope of this work. 

The changes occurring in the spectral radiance scattered by spherical 
grains of different composition and grain size distribution as a function 
of the phase angle (g) are calculated for homogeneous particles made of 
possible cometary material endmembers, including water ice, organic 
ice tholin, troilite, amorphous carbon and magnesium silicate. Fora 
given composition and phase angle, the single particle phase function 
p(g,A) modulates the intensity of the scattered spectral radiance. 

The theoretical trends of A,,,, aS a function of the scattering angle 
6=180° - gfor particles of different composition and three grain size 
distributions (radius 0.10-0.99 um, 1.0-9.9 pm and 10-50 pm) con- 
sidered in this study, are shown in Fig. 3. In the same figure, we show 
the distribution of /,,,, as measured from VIRTIS observations (Fig. 3f) 
selected during six MTP periods encompassing pre-perihelion (MTPO15 
and MTPO16), perihelion (MTPO20 and MTPO21) and post-perihelion 
(MTPO27 and MTPO28) orbital phases. These periods are chosen 
because they offer the widest spread in scattering angle necessary 
to discriminate among different compositions and grain size distri- 
butions. Our analysis shows that, with the exception of transparent 
grains, for example, water ice and organic ice tholin, the response of 
Opaque materials like amorphous carbon (Fig. 3a) and troilite (Fig. 3d) 
are characterized by aconstantA,,,,,at about 0.5 pm for the intermediate 
(1.0—9.9 um) and large (10-50 pm) grain size distributions across the 
60 < 8<150° scattering angle range explored by VIRTIS. These grains 
are therefore not matching the observed distribution of the A,,,,. Con- 
versely, opaque-material submicrometre grains show a significant shift 
of the A,,,,, towards the red range, especially for 9< 120°. 

Submicrometre grains made of magnesium silicate (Fig. 3b) show 
a linear increasing trend of A,,,,. from 0.56 um at 8 = 60° to 0.63 um at 
6=150°, whichis not compatible with VIRTIS data. In contrast, the A,,,,, 
of grains with radii larger than 1m is dispersed around 0.55 um, a value 
too red to reproduce VIRTIS data at > 100° and too blue for 8< 100°. 

Water-ice grains (Fig. 3e) offer the only solution matching VIRTIS 
data taken during the pre- (MTPO15 and MTPO16) and post-perihelion 
periods (MTP0O27 and MTPO28) at scattering angles 6 > 100°. During 
these phases, VIRTIS data are dispersed at A,,,, < 0.5 um, a behaviour 
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compatible with only the three water-ice grain-size distributions we 
have simulated. 

Apart water ice, organic material grains (Fig. 3c) are the second end- 
member showing spectral similarities compatible with 67P particles: 
the 1-10 pm and10-—50 pm grain-size responses are in fact characterized 
by alinear decrease of A,,,,, from @=60° to 120°, making them compat- 
ible with observations acquired during perihelion and post-perihelion 
periods. In contrast, submicrometre grains show an extreme red colour 
across the entire range of scattering angles: having A,,,, > 0.65 um, they 
are marginally compatible only with VIRTIS observations taken at low@ 
angles. If one assumes that during one month, for example, during the 
duration of one MTP period, the dominant composition and grain size 
is preserved, then the VIRTIS measurements taken during MTPO28 at 
6>130° could be compatible not only with water ice but also with large 
(10-50 pm) organic material grains. As shown in Fig. 3a, amorphous 
carbon submicrometre-size (O.1-0.9 um) grains are characterized 
by areddening of A,,,, for 8< 90°, similar to VIRTIS but with absolute 
values too low to exactly match the observations. Instead, the match 
is compatible for 90° < 8< 120°, making submicrometre amorphous 
carbon grains a good interloper between water-ice and organic grains. 
Finally, magnesium silicate (Fig. 3b) and troilite (Fig. 3d) submicrometre 
grains show A,,,,, trends at low scattering angles not compatible with 
VIRTIS data. 

Composition and grain size distribution of the dust particles were 
further analysed through the 0.5—0.8 pm spectral slope (Fig. 4) derived 
form VIRTIS data, which shows a peaked distribution at 8 = 100-110° 
(Fig. 4f) occurring in perihelion (MTPO20 and MTPO21) and post peri- 
helion (MTPO28) sequences. In our simulations, a similar behaviour is 
compatible with submicrometre-size grains made of organic matter 
(black points in Fig. 4e) or with greater than 1 um water-ice grains (blue 
and red points in Fig. 4d). The absolute value of the slope on the peak 
(about 3 um‘) measured at perihelion (MTPO20, cyan points) is similar 
to theoretical values for organic matter. Being transparent, water-ice 
particles are characterized by a peculiar behaviour: submicrometre 
grains show a blue colour response (negative slope, black points in 
Fig. 4d) not compatible with any VIRTIS observation, whereas large 
grains (blue and red points) are neutral to moderately red and could 
contribute to the VIRTIS observed slope at any scattering angle. Opaque 
grains with sizes greater than 1 um show an almost constant neutral to 
moderately red (0.5 ym") slope within the entire scattering angle range 
(see red and blue points for amorphous carbon in Fig. 4a and troilite 
in Fig. 4b). We remark that while the A,,,,, trends of magnesium silicate 
(Fig. 3b) grains are not compatible with VIRTIS data, their spectral 
slopes values (Fig. 4c) for radii greater than 1 um are remarkably similar 
for measurements taken during the post-perihelion period. 

Extended Data Table 3 shows the compatibility scheme of the A,,,, 
and 0.5—0.8 ppm spectral slope values for the different composition 
endmembers, grain radius and scattering angle compared with the 
VIRTIS data. 


Data availability 


The VIRTIS calibrated data are publicly available through the European 
Space Agency’s Planetary Science Archive website (https://archives. 
esac.esa.int/psa/). 
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Extended Data Fig. 1| Example of atypical VIRTIS-M observation of 67P 
nucleus and coma. Top left: visible colours RGB = (0.7, 0.55, 0.44 um) image, 
stretched to saturate the nucleus and enhance the visibility of jets inthe coma. 
Top right: tangent altitude image where the green mask corresponds tothe 
annulus containing all pixels witha tangent altitude between1land2.5kmfrom 
the limb. Bottom left: average radiance spectra as derived from official 
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pipeline (black curve) and after correction with Vega data (red curve). The 
fourth-degree polynomial fit to the corrected radiance is shown (cyan curve). 
The retrieved maximum emission wavelength on the fit is indicated by the 
magenta dashed line. Bottom right: corresponding //Fspectra normalized at 
0.5 um. The best-fitting slopes inthe 0.4—0.5 and 0.5-0.8 um ranges are 
indicated by blue and green dashed lines, respectively. 
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Extended Data Fig. 2| Vega star signal and derived VIRTIS responsivity. Top 
left: average spectrum of Vega, in DN, from observation V1_00402035638.QUB 
(red curve) and spectrum corrected for dark current, despiked and filter notch 
removed (green curve). Note the correlation of negative spikes in the Vega and 
average sky spectrum (black curve). Owing to the instrumental point spread 
function and spectral tilt, the signal is an average taken from 18 pixels. Top 
right: average spectrum of Vega derived from the four observations listed in 
Extended Data Table 2 after having applied a processing similar to the one 
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shown in the previous plot. The curve is averaged with a two-point running 
boxcar filter. The Vega flux** is shown in relative units (blue curve). Bottom left: 
VIRTIS responsivity derived from Vega signal averaged with atwo-point 
running boxcar filter (red curve) and nine points (blue curve) after 
normalization at 0.635 um above the standard responsivity value. Bottom 
right: comparison between standard pipeline (black curve) and Vega 
responsivities with a nine-point running boxcar filter (red curve). The ratio 
between thetwo responses is the blue curve. 
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Extended Data Fig. 3 | Rosetta spacecraft three-dimensional trajectory and 
solar phase angle variations with time. Top: Rosetta trajectory inthe 67P XYZ 
reference frame. Points along the X axis are shown starting from Rosetta’s 
position at 2015-01-131T23:28:53 (MTPO12) with an increment of 1km every 

20 min to improve visualization. The Yaxis is oriented towards the Sun and the 
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Zaxis is perpendicular to the orbital plane. The red line indicates the position 
of the nucleus along the X axis. The integrated radiance as measured on each 
observation is reported according to the colour scale. Bottom: variation of the 
solar phase angle (Sun—nucleus centre—Rosetta) during the mission. 
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Extended Data Table 1| Rosetta’s calendar MTP periods dates and heliocentric distance 


Period Heliocentric distance (AU) 





Extended Data Table 2 | List of VIRTIS observations of the Vega star 


Observation Start Time End Time (Band, Int. 
Sample, Time (s) 
Line) 
1_00402035638.QUB | 2015-09- 2015-09-28T | (432,256,29 50 
28T04:35:17.781 | 05:04:11.405 


V1_ . : : ( ) 

V1_00403369562.QUB | 2015-10- 
19T04:03:18.373 | 04:34:12.013 
131T15:36:39.065 | 16:03:32.757 





Article 


Extended Data Table 3 | Summary of spectral indicator compatibility 
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10-50 
0.10-0.99 
Troilite 1.0-9.9 
10-50 
0.10-0.99 
Water Ice 1.0-9.9 
10-50 
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The table shows the compatibility of spectral indicators (A,,,,, O.5-0.8 um spectral slopes) between VIRTIS observations and Mie scattering simulations for a given composition, grain radius 
range and scattering angle. The compatibility between observations and scattering simulations is colour coded according to the key shown on the top right. Each cell is divided in four fields 
showing J,,,, and to 0.5-0.8 um slope values during pre-perihelion, perihelion and post-perihelion epochs (bottom right key). 
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